Structural Relaxation of a Gel Modeled by Three Body Interactions 
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We report a molecular dynamics simulation study of a model gel whose interaction potential is 
obtained by modifying the three body Stillinger- Weber model potential for silicon. The modification 
reduces the average coordination number, and suppresses the liquid-gas phase coexistence curve. 
The low density, low temperature equilibrium gel that can thus form exhibits interesting dynamical 
behavior, including compressed exponential relaxation of density correlations. We show that motion 
responsible for such relaxation has ballistic character, and arises from the motion of chain segments 
in the gel without the restructuring of the gel network. 



Gels are low density disordered networks of interact- 
ing molecules that are structurally arrested and capable 
of sustaining weak stresses. They are ubiquitous in na- 
ture and among man-made materials and are composed 
of a diverse range of materials such as polymers, silica, 
or colloidal particles. Depending on the life time of the 
bonds between the basic units of the network, they can ei- 
ther be chemical gels, or reversible, physical gels, the lat- 
ter displaying complex dynamics. In particular, colloidal 
gel formers exhibit intricate dynamic behavior in equi- 
librium, as well as in noncquilibrium aging conditions, 
can form arrested states, and have been the subject of 
a considerable number of experimental, theoretical and 
simulation studies J^, [1J|, i,B H 0, 0, S EI El EI El, 
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One reason for the interest in colloidal gels is that 
these systems permit access to the glassy state via several 
mechanisms: Cluster aggregation [3, [l9(, structural ar- 
rest in the dense phase following phase separation [I EH , 
or crossing of a glass transition line from an equilibrium 
fluid to an arrested state d 15 1. For the occurrence 



of this latter scenario it is necessary that upon cooling 
the system does not enter the liquid-gas coexistence re- 
gion E EI EI EI El, i.e. one seeks systems for which 
the coexistence region is at low temperatures, T, and 
densities, p. One possibility to achieve this is to choose 
a "maximum valency" interaction, in which each parti- 
cle can interact only with a restricted (small) number of 
particles [I EI EH ■ In the following we will show that 
a very simple model involving three body interactions is 
also able to generate a coexistence region that is located 
at low T and p and which therefore allows one to probe 
easily the interplay of phase transformations and dynam- 
ics in molecular dynamics simulations in such systems. 

A further intriguing property of colloidal gels is the fact 
that their relaxation dynamics can be compressed 0, 01 j 
i.e. the time correlation functions decay faster than an 
exponential, in stark contrast to structural glasses at 
higher densities for which one usually finds a stretched 
exponential relaxation. The microscopic origin of this 



fast relaxation is not well understood, and various mech- 
anisms have been proposed to explain it d EI H3] ■ We 
present analysis that shows that for our model, com- 
pressed exponential relaxation arises from the ballistic 
motion of chain segments in the gel without restructur- 
ing the gel network. 

The model we consider is a modification of the po- 
tential proposed by Stillinger and Weber (SW) for the 
description of silicon j^. Particles interact via a sum 
of two and three body interaction terms, v = i'2(r) + 
A«3 (r, 9) [23| , where r denotes interparticle distances and 
6 the angle formed by three particles. A determines the 
strength of the three body interaction which depends on 



the angle via a term proportional to (cos / 



with 1 



determining the most preferred angle. Thus by varying A 



24 1 and a we can tune the locally preferred arrangement 



of the particles. 

We have performed constant temperature, volume 
molecular dynamics (MD) simulations (using a constraint 
that conserves kinetic energy) with 4000 particles, using 
the method proposed in [H 0, [23] to efficiently com- 
pute three body interactions. Gibbs-Ensemble-Monte- 
Carlo (GEMC) simulations [28| are performed to obtain 
liquid-gas coexistence curves have been performed with 
2000 particles. All results are reported in reduced units 
for the Stillinger- Weber potential [22 1. 

Figure Q] shows the coexistence curves obtained for var- 
ious combinations of A and a. We see that with increas- 
ing A or a, liquid-gas phase coexistence gets shifted to 
smaller temperature and density ranges, analogous to 
the observations in [I EI EH ■ In the following we will 
fix A = 10.0 and a = 1.49. For this choice the struc- 
ture of the system at low T and p is given by quasi-one- 
dimensional chains of particles, interconnected by three 
coordinated junctions. At low T bond breaking becomes 
extremely difficult and consequently a reliable estimate 
of the coexistence curve via GEMC is no longer possi- 
ble. Based on MD runs where we observe signatures of 
phase separation, we indicate in Fig. Q]the region where 
we expect phase separation (shaded area) . Also included 
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FIG. 1: Phase coexistence curves calculated from GEMC sim- 
ulations for (a) A=21, a=l/3, (b) A=25, a=l/3, (c) A=25, 
a=l/2, (d) A=10, a=1.00 and (e) A=10, a=1.10. With the 
increase of value A or a the phase coexistence curve is sup- 
pressed. The shaded area indicates the expected coexistence 
region for A = 10.0 and a = 1.49. Also shown is the percola- 
tion line for these values of A and a. 



in the graph is the percolation line which indicates the 
density and temperature range (to the lower right of the 
percolation line) where we may expect gel-like structural 
arrested states. (The bend in the percolation line at low 
temperature is due to phase separation [27j]-) In the fol- 
lowing we will study the relaxation dynamics of the sys- 
tem for p = 0.06 and from Fig. |T] it is clear that at this 
density phase separation will not play a role. 
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FIG. 2: The collective intermediate scattering function 
F(k, t) (left panel) and self intermediate scattering function 
F s (k,t) (right panel) at T = 0.03, for a range of wave- vectors 
k. 

In order to characterize the relaxation dynamics of the 
system we consider the normalized collective, and self 
intermediate scattering functions, defined as F(k, T) = 
S^k)- 1 £j,i(expHk- ( Vj (t) -r;(0))]}, where S(k) is the 
static structure factor, and F s (k,T) = N^ 1 £ (exp[— «k- 



(fj(t) — r,-(0))]}, respectively. The time dependence of 
these correlators is shown, in semi-log plots, in Fig. [2l 
for the low temperature T = 0.03, i.e. well below the 
percolation line, and various values of the wave-vector 
k. At this T the relaxation dynamics of the system is al- 
ready very sluggish and hence we deal here indeed with a 
glass- forming system (see Fig. [4] in which one observes a 
strong ch ang e of the relaxation times with temperature, 
and Rcf. [27f for a detailed discussion). Since we look for 
compressed exponentials, we plot the data as a function 
of t /t and t /t s , where r and r s are the relaxation times 
defined by requiring that the correlator has decayed to 
e _1 of its initial value. We see that the two time cor- 
relation functions display remarkably different behavior: 
For intermediate wave-vectors F(k,t) curves downward 
and can be fitted well by the Kohlrausch- Williams- Watts 
(KWW) function A exp(-(i/V)' 3 ) with (3 w 3/2, the so- 
called compressed exponential (CE). [If k is very small 
or very large the decay is even faster, i.e. > 1.5 (not 
shown)]. Such a behavior (/3 w 3/2) has been observed in 
experiments of slowly relaxing gels [1, [f| , and analyzed 
theoretically using a stress relaxation model [2(| . On the 
other hand F s (k,t) shows, like most other glass-forming 
systems, a stretched exponential, i.e. (3 < 1.0 for small 
and intermediate wave-vector, and a compressed expo- 
nential at large wave-vectors. Such a behavior has also 
been observed in [l6j], and interpreted as arising from 
the averaging of ballistic motion of particles which form 
a part of chain segments of varying lengths in the disor- 
dered percolating network. 

To analyze further the nature of the relaxation, we 
determined the KWW-exponent (3 by fitting the correla- 
tors to a KWW-function. Since the time correlators ex- 
hibit different regimes of decay, it is necessary to choose 
a meaningful and consistent procedure for obtaining /3 
and we choose to fit the curve in the time window in 
which the correlator is between 0.9 and 0.1. This choice 




FIG. 3: Wave-vector dependence of the KWW exponent f3 
for F(k,t) and /3 3 for F s (k,t) for four different temperatures. 
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avoids the (trivial) ballistic regime at very short times 
and instead focuses on the relaxation regime seen at in- 
termediate times. The wave-vector dependence of the 
so obtained KWW-exponents is shown in Fig. [3] for four 
different temperatures. The (3 values for F(k,t), shown 
in the top panel, are seen to be always 2.0 for T = 5.0, 
indicating that ballistic motion dominates the decay (as 
discussed in the context of gels in at all wavelengths. 
While this is expected for large k, we note that the be- 
havior at small k is a result of the low density of our 
system which leads to significant decay of collective den- 
sity fluctuations even on large wavelengths through non- 
diffusive motion of particles. If T is decreased (3 shows 
a minimum at intermediate values of k and the width 
of this minimum broadens with decreasing T, suggesting 
that on intermediate length scales the decay mechanism 
is distinct. The typical values of (3 in this minimum are 
around 1.3 — 1.6, i.e. similar to the values that have been 
found in the experimental systems or in the theoretical 
calculations. 

The p s values for F s (k,t), shown in the bottom panel, 
are found, for the highest temperature T = 5.0, to change 
from 2.0 at large k towards 1.0 at small k, which we inter- 
pret as the expected crossover from ballistic to diffusive 
decay. Again, for low temperatures, we find superposed 
on this overall trend an intermediate regime, in which 
the dynamics becomes "stretched", i.e., the decay of the 
self motion becomes slower than exponential. Thus from 
this figure we can conclude that the self and collective 
density correlation functions exhibit complex behavior, 
that is non-trivial, and different from dense fluids. 
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FIG. 4: The wave-vector dependence of relaxation time from 
the area of F(k,t) and F 3 (k,t) for temperature 0.028, 0.06 , 
0.30 and 5.00 at density 0.06. The filled symbols are for t 3 
and opaque symbols indicates r. 

The relaxation times, r A (fc) and T A (fc), obtained from 
calculating the area under F(k, t) and F s (k, t), are shown 
in Fig. H] as a function of k for different temperatures. 
Since for ballistic motion one expects the relaxation time 
to be proportional to A: -1 , we show r A (k) and T A (fc) mul- 



tiplied by k. From the figure we recognize that at high 
and intermediate T this scaling gives indeed horizontal 
line, showing that the motion can be interpreted as bal- 
listic. Furthermore we see that the self and collective 
relaxation times track each other for all k. At low T 
and small k the curves are no longer horizontal, indicat- 
ing that there is a significant non-ballistic component . As 
we shall see below, at these low temperatures, F(k, t) has 
a significant long time relaxation that is clearly distin- 
guishable from an intermediate relaxation process which 
we shall identify with CE behavior. Furthermore there 
is a strong decoupling at low k in that t a exceeds r A by 
a large factor. 
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FIG. 5: (a) The F(k,t) from molecular dynamics (MD), and 
Langevin dynamics (LD) for different damping constants 7. 
Also shown is F(k, i) from constrained MD simulation de- 
scribed in the text, and the KWW fit to the MD curve for 
intermediate times, (b) rk vs wave-vector, k for a range of 
temperatures, showing that at low k, t ~ 1/k. 



We now investigate the nature of the compressed ex- 
ponential relaxation of F(k, t) (shown for T = 0.04, k = 
0.15 in Figure 5(a)), which we observe on intermediate 
time scales. At low temperatures, CE relaxation accounts 
for a substantial part of the decay of F(k, t) for a wide 
range of intermediate k values. Since at these low tem- 
peratures the average life time of the bonds is longer 
than the decay time of F(k,t) in the CE regime 27| . 
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we expect the CE to be associated with the floppy dy- 
namics of chain segments in the transient gel network, 
without network restructuring playing any role. In or- 
der to analyze the motions that are relevant, we there- 
fore compare the F(k, T) obtained in MD simulations to: 
(a) MD simulations with the imposition of a constraint 
that prevents bond breaking and formation of new bonds. 
This is accomplished by identifying bonded neighbors in 
the initial configuration we consider, and adding a suit- 
ably parametrized barrier potential of gaussian form to 
the two body part of the S-W potential, (b) Langevin 
dynamics simulations, to study the role of microscopic 
dynamics. For the Langevin dynamics we have used a 
predictor-corrector integrator [29(, and the damping co- 
efficient is tuned to span the range from very small damp- 
ing, 7 = 0.004 (corresponding to MD) to strong damping, 
7 = 0.4. In Fig. [5] (a) we show F(k, t) from these differ- 
ent simulations. 

Comparing MD with the constrained MD results, we 
see that the regime of CE dynamics is essentially unal- 
tered by the imposition of the constraint not to break 
or form bonds. This shows clearly that the CE dynam- 
ics arises from the dynamics of the non-restructuring gel 
network. However, at longer times the relaxation dynam- 
ics of the constrained MD is essentially frozen, indicating 
that long time relaxation in the MD, cleanly separated 
from the compressed exponential decay, arises from net- 
work restructuring, a result that is also confirmed by the 
time dependence of the mean squared displacement of 
the particles (not shown) [27j |. 

In Fig. [5] (b) we show r x k where r(fc) is obtained by 
fitting F(k,t) by a KWW function in the intermediate 
time window displaying CE relaxation, as shown in Fig- 
ure 5 (a) (note that these t values are different from those 
shown in Fig. [4] which are obtained from the area under 
F(k, t)/S(k)). The near constant value of r x k for small 
k corroborates the ballistic origin of the compressed ex- 
ponential relaxation, consistent with predictions [3. |20|. 

For the Langevin dynamics we see that for small and 
intermediate values of the damping coefficient 7, the cor- 
relator tracks the one from the MD and thus a CE will 
be observed. However, for large damping the shape of 
the curve is very different from the one of the MD and 
no CE is seen anymore. Thus we see that the dissipativc 
dynamics, relevant for example for real colloidal gels, will 
not show a CE dynamics. Therefore we can conclude that 
the CE seen in those systems is likely due to the aging 
dynamics. 

In conclusion, we have proposed a model system which 
allows the simulation and study of gel forming fluids un- 
der equilibrium conditions, by suppressing the liquid-gas 
phase coexistence curve to an arbitrarily small tempera- 
ture and density window. At low densities and temper- 
atures the structural and dynamical features show many 
similarities to the one of experimental systems. In par- 
ticular wc find an intricate behavior of the density cor- 



relation functions, including compressed exponential re- 
laxation of the collective intermediate scattering function 
with a compressing exponent that depends on tempera- 
ture and wave- vector considered. The motion responsible 
for the compressed relaxation is found to have ballistic 
character and to arise due to the motion of chain seg- 
ments in the gel without the restructuring of the gel net- 
work. 
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